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We consider a general class of discrete unitary dynamical models on the lattice. We show 
that generically such models give rise to a wavefunction satisfying a Schrodinger equation in 
' the continuum limit, in any number of dimensions. There is a simple mathematical relation- 

, ship between the mass of the Schrodinger particle and the eigenvalues of a unitary matrix 

' describing the local evolution of the model. Second quantized versions of these unitary 

models can be defined, describing in the continuum limit the evolution of a nonrelativistic 
quantum many-body theory. An arbitrary potential is easily incorporated into these sys- 
tems. The models we describe fall in the class of quantum lattice gas automata, and can 
be implemented on a quantum computer with a speedup exponential in the number of par- 
ticles in the system. This gives an efficient algorithm for simulating general nonrelativistic 
' interacting quantum many-body systems on a quantum computer. 



I. INTRODUCTION 

There are many situations in physics where a continuous system obeying a particular set of equations at a 
macroscopic scale can be modeled by a discrete microscopic system obeying a very simple set of local rules. 
For example, in equilibrium statistical mechanics, simple lattice models such as the Ising model capture the 
behavior of generic clasaesr-pf critical systems at large scales. Another interesting class of discrete systems 
are lattice gas automataErEl; these models describe systems of particles moving about on a lattice, obeying 
simple collision rules which conserve quantities such as mass and momentum. In the macroscopic limit, these 
systems obey Navier-Stokes or other hydrodynamic equations of interest. 

In the quantum domain, there are also examples of discrete microscopic systems which capture interesting 
macroscopic behavior. Lattice-gauge theories (see for example Creutzo) give an approach to studying the 
partition function and spectra of quantum field theories by mapping these theories to statistical mechanical 
ensembles. There are, however, few discrete models for describing the dynamical evolution of quantum 
systems which preserve important features such as unitarity. An example of a quantum system for which a 
unitary discrete model is known is the, Dirac equation describing a relativistic particle moving in one spatial 
dimension. As shown by FeynmanBla, this system can be described by a simple microscopic model of a 
particle moving on a ID lattice according to a simple local rule which essentially corresponds to a unitary 
form of random walk. A straightforward attempt to realize a Dirac equation in more than one spatial 
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dimension as a form of unitary random walk cannot succeed. By using operator splitting methods, however, 
it was shown by Succi and BenziQ that a sequence of random moves along single axes, alternating with 
transformations which diagonalize each of the Dirac matrices in turn, can give an analogous construction 
in higher dimenaiai|is. The discrete model for the (1 + l)-dimcnsional Dirac equation has been of renewed 
interest recentlyOLl, due partly to the possibility of simulating such, unitary microscopic discrete systems 
by quantum computers. In particular, it has recently been suggestecB that a simple quantum lattice model 
can be constructed which describes the motion of a system of many particles moving according to the one 
dimensional Dirac equation. 

In this paper we consider a class of models closely related to the ID Dirac lattice model, which give rise to 
a nonrelativistic single-particle Schrodinger equation in an arbitrary number of dimensions. In these models, 
the time development rule is given by a single local, unitary transformation matrix. Thus, we are essentially 
considering the motion of a single particle under a unitary random walk process. For this class of models 
we show that the macroscopic equation of motion satisfied by the wavefunction corresponding to particle 
density is the Schrodinger equation. We show that such nonrelativistic models can be constructed for an 
arbitrary number of spatial dimensions. We also show that an arbitrary potential can easily be included in 
these models. 

It is natural to generalize from the single particle models to a second quantized many-body system. Such 
a model could be implemented very efficiently on a quantum computer, so that the number of computational 
steps necessary to simulate a single time step would only depend upon the size of the lattice and would 
not depend upon the number of particles in the system being simulated. Thus, our results could be used 
to efficiently simulate an arbitrary nonrelativistic interacting quantum many-body system on a quantum 
computer exponentially faster than the same calculation could be performed on a classical computer. Such 
a system would be an ideal example of quantum computing, since the computing elements could be built 
from a system of spin-1/2 particles on a lattice obeying simple local unitary time evolution rules. 

The principal .-difference between our models and the ID Dirac lattice model (and its generalization by 
Succi and BenziQ) is that in the Dirac model, the unitary evolution rule satisfied by the wave function or 
particle at each time step is infinitesimally close to the identity transformation. As the lattice spacing e goes 
to 0, the unitary transition matrix is of the form S = 1 + ieM, where M is Hermitian. In our models, we take 
the transition matrix to be independent of the lattice spacing. This form of a time development equation 
makes the system nonrelativistic, but allows for a formulation in an Jixbitrary number of dimensions. A 
closely related model was considered recently in one spatial dimensionE2l, where simulations were shown to 
be consistent with emergent behavior corresponding to a Schrodinger equation. In this paper we prove that 
this is the general behavior of such models, giving a simple algebraic relation between the transition matrix 
and the mass of the nonrelativistic particle. We develop such models for the Schrodinger equation in an 
arbitrary number of dimensions. 

In the first part of this paper we will consijder lattice models for single-particle motion. These models are 
essentially unitary lattice-Boltzmann modelstil. In the latter part of the paper we generalize to many-body 
systems and discuss how a lattice of quantum computing elements could be used to describe the motion 
of a large number of nonrelativistic quantum particles. We conclude with a simple numerical check of the 
analytic description of a sample model. 

As a simple example of the type of system considered in this paper, consider the lattice-Boltzmann model 
with a configuration space defined by two complex fields, i^iix, t) and ip2{x, t), taking independent values on 
a lattice with one spatial dimension x and one temporal dimension t. Define the dynamics of this model to 
obey the equations 

Mx + l,t) = ^[{1 ~ i)2/ji{x,t - 1) ~ {1 + i)Mx,t ~ 1)] 



- 1, i) = - [(1 - i)^2{x, t - 1) - (1 + i)M^, ^ - 1)] • 

These equations give a unitary time evolution to ijj. To understand how ip evolves in a continuum limit, we 
can expand the equations of motion through 4 time steps, giving for example 
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V'i(a;,i + 4) = ^ [--iPiix ~ 4:,t) + 3iJi{x - 2, t) + 'iPi{x,t) + tPiix + 2, t)] 

+ ^ [V'2(x - 2, t) - t) -Mx + 2, t) + Mx + 4, i)] • 

Taking a continuous limit as the lattice spacing scales as e in the x direction and in the t direction, we 
find the diS^erential equation 

dtMt) = 'flMt) 

A similar equation holds for -02, and so it follows that 

Thus, we see that the total amplitude '^{x,t) — tpi{x,t) + tp2{x,t) satisfies a Schrodinger equation. As we 
shall demonstrate, this is the generic behavior of a unitary Boltzmann model with a fixed time development 
rule. 

We introduce the Schrodinger model in Sec. || by presenting the one-dimensional case. The model is 
generalized to Cartesian lattices of arbitrary dimension in Sec. in this section we also discuss the inclusion 
of a potential. In Sec. [V, we discuss how the one-particle models can be generalized to construct a quantum 
lattice-gas model of many nonrelativistic particles. In Sec. we give the results of a simulation of a single 
free nonrelativistic particle in 2D, comparing numerical results with the theoretical framework presented 
here. The appendices give explicit formulas for models in 2D and 3D on a Cartesian lattice. 



II. SCHRODINGER EQUATION IN ONE DIMENSION 



In this section, we consider unitary lattice-Boltzmann models describing the evolution of a single particle 
in one dimension. Keeping the collision operator fixed in the scaling limit, we show that a very general class 
of microscopic models give rise in the continuum limit to a Schrodinger equation. 

We define the model on a lattice given by points x — en where n is an integer. The lattice can be taken to 
either have periodic boundary conditions or to be of infinite extent. The state of the system at a fixed value 
of the time parameter t is described by a wave function 4'k{x, t) which depends upon the discrete position x 
and an "internal" index k taking values from 1 to to, labeling possible particle velocities at the lattice site x. 
As in lattice-Boltzmann models, at each discrete time step the various components of the field at each site 
undergo a local unitary "collision," and then the jth component of ip{x,t) propagates along the jth lattice 
vector Cj to the new site x + Cj to yield the new state of the system at time t + At. We consider only linear 
processes, so this interaction can be specified by an TO-by-m scattering matrix S. 

We take the continuum limit of the theory by scaling e ^ 0, where At ^ e^. In this limit we will find 
that the discrete equation describing the dynamics of ip becomes a continuous differential equation, which 
we identify as the Schrodinger equation. 

In this section we will assume that each lattice site has two associated possible particle velocities, cor- 
responding to right- and left- moving particles. We will also assume that the dynamics is symmetric under 
right-left refiection. More general models can easily be analyzed using a similar formalism. 

The equation of motion for the model reads 

ipk{x + eck.t) = Skjtpj{x,t- At) 

where fc = 1,2 correspond to right and left moving particles, so that ci — +1,C2 — —1. The quantum wave 
function ipk{x,t) is normalized so that 

J2\Mx,t)\'^i (ii.i) 

x,k 
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for all t. The matrix Skj is a two- by-two matrix which is unitary so as to preserve the condition Eq. ( ]II.1| ). 
We begin our analysis by transforming to the wavefunction 

tpix, t) — S'^0(a;, t), 

where t = t/At. We can then expand the difference (j){t) — ip{t — 1) in the infinitesimal parameter e, to get 
(p{t) - (j){t - 1) = -eS-''CS^d^(j) - —S-^C'^S^dl<p 0{e^) 



where C is the two-by-two matrix 

C 



I 

-f 



Because we are assuming that the interaction described by the matrix S is invariant under reflection, S 
must be of the form 



S = 



a b 
b a 



where a and b are complex numbers. Because of unitarity we have |ap -I- |6p = 1.5' can be put in diagonal 
form by writing 

S = X-^DX 

where 

W 1 1 



V2V1 -1 



We can redefine S and up to a phase, so that without loss of generality we can take 

D = 



pi 
1 



where is a complex number with magnitude 1 (/i/i* = 1). We then have 

2 ^ M - 1 M + 1 
If we write X(j) — r;, then we have 

77(i) - r^it - 1) = -eD-'{XCX-^)D''d.^ri - —D-^XC^X-^)D^dl7^ + ©(e^) 

At this point we would like to scale the time step as a power of e so that this equation can be written 
as a differential equation in time. However, there is a difhculty which arises due to the fact that there are 
two relevant time scales involved in the dynamics of r/. There is an order-e change to rj at every time step; 
however, this order-e term has a phase angle which rotates at every time step. Thus, the order-e dynamics 
average out after a large number of time steps, so that the time-averaged rate of change of rj actually goes 
as e^. The dynamics we are interested in are independent of the short-term order-e fluctuations, so we must 
perform another transformation to remove these effects. With this goal in mind, we write 

= at) + epit) 

where 
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p{t) - p{t - 1) = -D-^XCX-^)D^d^(:- 

This equation is solved by 

p{t) = D--GD-d,C, 

where 

G - DGD-^ = -B = -{XCX-^). 

This can be solved for G as long as the only nonzero entries Bij appear where the i and j eigenvalues of D 
are different. We can now write a final dynamical equation for C,. 

C(i) - C{t - 1) = -e'^D-^BGD^dlC, - -D-''{XC'^X-'^)D''dlC + 0{e^). 
If we assume that the unit of time scales as e^, we have the continuous dynamical equation 

dtC = -D-'BGD-'dlC - ^D-^{XC^X-^)D''dlC- 



We can now substitute the known matrices X, C, D to compute 

B = XGX-^ -- 



1 

1 



xc'^x-'^ = 



1 
1 



Using these matrices we have 



1 ' ' 



BG + -XC^X-^ = ( 2 ^i-/'* ^ ^ 



2 



Writing p = cos 9 + i sin 9, we have 

1 1 — cos ( 



(l-cos^)2 + sin2 6l 2 2(1 -cose*)' 
Thus, the dynamical equation for ( becomes 

9tC = * f 4 _D die 

where 

m = cot 9 — CSC ^. 

The equation for the first component of C is thus precisely a Schrodingcr equation for a particle moving in 
one dimension with mass m. To leading order, ({t) is related to tp through the sequence of transformations 
described above, so that 
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(it) = D-^Xm + 0{e) 
The first component of C,[t) is tlierefore given by 

Vl/ = ClW = ^(V^lW+^2W), 

and this satisfies the Schrodinger equation in the continuum limit, 

1 



2m 



Note that by taking fi = —i we get m = 1, giving precisely the example discussed in Sec. ||. We shall 
demonstrate in Sec. Ill that, in an analogous fashion, in higher-dimensional theories the sum of the wave 
function components forms a scalar quantity which satisfies a Schrodinger equation. 



III. SCHRODINGER EQUATION IN DIMENSIONS D > 1 



In this section we derive the general form for the continuum limit of the dynamics for a unitary lattice- 
Boltzmann model with fixed collision matrix on a lattice with any number of dimensions. Specializing to the 
case where the lattice is Cartesian and the collision rule is invariant under discrete rotations, we find that a 
generic collision rule gives a Schrodinger equation in any dimension d. 



A. General form of dynamical equation 



The analysis of the continuous equations of motion in d dimensions proceeds in a fashion very similar to 
the discussion in the previous section. We assume that the lattice contains a set of points x, and that at 
each lattice site there are particle velocities labeled by k, corresponding to velocity vectors in the lattice. 
Denoting spatial indices by a, we denote the ath component of the fcth velocity vector by c^. The dynamics 
of the lattice-Boltzmann model are described by the equation of motion 

ipki'x. + eck,t) = S'fejVj(x,i - 1) 
where S is unitary. Transforming as before 

we have 

(f>{t) - (P{t - 1) = -eS-^C"5^5a0 - —S^^CC'^S^da.d^cl) + 0{e^) 
where the diagonal matrices C" are given by 

= diag(c?,...,<), 

with c" being the ath spatial component of the jth lattice vector. Writing S — X^^DX, X4> — f] we have 
r]{t) - 7]{t - 1) = ~eD-''{XC°'X-^)D''d„r] - -D'^ {XCCf^ X-^)D^d^dpT] + 0{e^) 

We write 



Vit) = at) + ep{t) 
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where 

p{t) - p{t - 1) = -D-'{XC"X-^)D^daC 

This is solved, as before, by 

p{t) = D--G"D^dc,C, 

where 

G - DGD-^ = -{XCX-^) = -B 

Again, this can be solved for G as long as the only nonzero entries Bij appear where the i and j eigenvalues 
of D are different. The resulting continuum equation for jj is 

dtT] = -D-''B°'G'^D''dad0r] - ^D-"'iXG°'C^X-'^)D''dad/3r] 
This is the general form of the dynamical equation for a unitary lattice-Boltzmann model. 



B. Schrodinger equation in d dimensions 



We now specialize to the case where the lattice is Cartesian, so that there are 2d possible particle velocities 
at each lattice site, corresponding to vectors of magnitude +e, — e in each of the d directions. We choose the 
collision matrix S to be invariant under the symmetry group of the lattice. Wc will show that gcncrically the 
continuum limit of the equation of motion is a Schrodinger equation, just as we found for a general collision 
matrix in ID on the Cartesian lattice. 

The constraint that 5* is invariant under discrete rotations and reflections is actually quite a strong condi- 
tion. The 2(i-dimensional space of velocity vectors transforms under a linear representation of this discrete 
group. This representation contains only 3 irreducible representations, which allows us to determine D up 
to 3 distinct eigenvalues. Because of the symmetry constraint, we can always diagonalize S by the matrix 



/ 



X 



V2d. 





1 



1 

V2 





1 _ 1 

X J. 

2V3 2V3 



1 









2 

2\/3 



V — !— — !— 

\ 2V3J 2v^ 2V5J 



1 

V2d 






1 






1-d 



1 




1 


!_ 

V2 





1 _1 

X X 

2\/3 2\/3 



1 









2 

' 2V3 



1 






1 

V2d 






-73 






2Vd2 2VS2 2v^ 2^32 2v^ 2V5J / 

where ^2 = d{d — l)/2. The rows of this matrix consist of the 3 groups of vectors in the irreducible 
representations of the rotation group mentioned above. The first row is the normalized vector (1, 1, ... , 1). 
The next d rows are normalized versions of the vectors c" with +1 in position i and —1 in position i + d. 
The last d — 1 rows are vectors with equal components i and i + d, subject to the condition that the sum of 
the components vanishes. This matrix puts S in the diagonal form 



v 












D = XSX-'^ = 








V 
A 



\0 ••• ••• A/ 
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where the eigenvalue v appears d times and the eigenvalue A appears d — 1 times. By a simple phase 
redefinition we can choose v = 1. We will furthermore set A = — 1, which as we shall see will give rise to a 
Schrodinger equation in the continuum limit. As we shall discuss later, however, any value of A ^ gives a 
Schrodinger equation; we use the A = — 1 condition merely to simplify the presentation. 

With the stated conditions on D, we can compute S. We find that all elements of S are equal to i^j^, 
except the matrix elements Si^i+n and Sj+d^j, which are equal to — 1. Thus, 

^ij - °0,\i-j\-d 

At a microscopic level, this collision matrix gives an equal amplitude for a particle to move in every direction 
other than directly backwards. To check the unitarity condition, we verify 

and 

^ ^ Ad^ 4rf2 Ad^ 

We can now proceed to calculate the other matrices needed for the dynamics. There are d matrices B°'. 
For a particular value of a, we find that all matrix entries vanish except those in the {a + l)th row and the 
(a + l)th column, which are given by 

{B'^)a+i,i = {B'^)w = (t^' 0'^""'' ^(")' + 1), . . . , r(rf - 1)) , 

where by 0^° we denote a sequence of k O's, and where we have defined 

r{a) = 



yja{a + 1) 

We can now immediately compute G", which has nonzero elements in the same positions, given by 



y/a-l 1 1 , 1 



We are now interested in computing the differential equation describing the continuum limit of the dynamics 
of the first component of Q, which we will denote hy ^ . As before, we define 

*(x,f) = Ci = ^ (^E^^(^'*)) +^(^)- 

To compute the dynamics of we need to know only the first rows of the matrices B'^G^ and XC"C^X"^. 
Prom the above expressions, we find that the first row of B"G^ is given by 

= (-;^- -r<«»' - • -'i 

From the above form of X, we see that the first row of XC°'C^X~^ is given by 
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(XC"C'^X-i)o. - S^f"^ (^_L, o'^+"-2, r(a), r(a + l),...,r{d- 1) 

Thus, the first row of the combined matrix is 



where 



m = d(cot 6 — CSC 9) 



(III.2) 



ith 



cost 



I sm( 



(IIL3) 



As a resuit, we obtain the differential equation describing the dynamical evolution of ^> in the continuum 
limit, 



1 



(IIL4) 



which we recognize as Schrodinger evolution in d dimensions. In Appendices ^ and ^ we work out the 
specific cases of ci = 2 and d = 3 in detail. 

Note that had we chosen the eigenvalue A differently, the difference would have appeared in the last d — 1 
rows and columns of G. Following through the computa tion, we find that the only change would be that 
new terms would appear on the right hand side of Eq. ( 111.4 ), proportional to the eigenvectors of S with 
eigenvalue A. However, these terms would have had a phase (/iA*)"^, and thus would have averaged out in 
the continuum limit as long as /i 7^ A, making no change to the final result Eq. ( 111.4 ). Thus, we reach the 
conclusion that for any collision rule invariant under the lattice symmetry group, so long as fi is distinct 
from the other eigenvalues of the collision matrix, the resulting continuum dynamics for the total amplitude 
^' are governed by a Schrodinger equation. 



C. Inclusion of a potential 



In general, we can easily include an arbitrary potential V^(x) by including a position dependent phase in 
the transition matrix S. If we perform the above analysis for a model with transition matrix 

^(x) = exp(-ieV(x))5 

where 5 is a spatially invariant matrix such as discussed above, then the general form of the dynamical 
equation becomes 

dtv{^,t) = -D--B^Gf'D-d^dpv{^,t) - ^D--{XC"C''X-')D^d^d0f]{^,t)-tVi^)r]{x,t). 
This becomes for the total amplitude ^, the Schrodinger equation in an external potential 

dt^i^, t) = V dl-^[^, t) - zy(x, i)vi/(x, t) 



9 



IV. MANY PARTICLES: QUANTUM LATTICE GAS AUTOMATA 



We now consider models in which multiple particles move independently according to the Schrodinger 
equation in d dimensions. One way of simulating the motion of n particles in d dimensions is to introduce 
extra degrees of freedom for each particle. Thus, for example, we could model the motion of 2 particles in 
one dimension by the lattice-Boltzmann model 

ipikix + eci,y + eck,t) = SiiSkjipijix,y,t - 1) (IV.5) 

where x,y are the positions of the two particles, i,k are the internal indices specifying their directions, and 
5 is a 2-by-2 matrix for unitary Schrodinger evolution in one dimension, as discussed in Sec. ||. Notice that 
this dynamics is equivalent to that of a single particle moving in two dimensions. 

In a similar fashion we can describe models where n particles move in d dimensions, by constructing 
a unitary lattice-Boltzmann model in nd dimensions. It is straightforward to incorporate an arbitrary 
interparticle potential in this formulation; the potential is a function of the particle positions and can be 
included as discussed in Section III.C. 

This gives a procedure for simulating an interacting nonrelativistic quantum many-body system on a 
classical computer. Although this may give a useful algorithm for systems containing only a few particles, 
if we wish to simulate the motion of a large number of particles using the method just described it is clear 
that the number of calculations needed to perform even one time step of the evolution become rapidly 
intractable. For example, simulating the motion of 20 particles in three dimensions on a lattice of side length 
100 would take on the order of 10^^" calculations per time step, beyond the capacity of any imaginable 
classical computer. j-, 

However, the technology of quanfMm computingi3 presents a paradigm in which such calculations can be 
done. We will now describe a way in which the above algorithm can be implemented on a quantum computer 
with a speedup exponential in the number of particles. In fact, it is natural to simultaneously perform the 
calculation for all numbers of particles which will fit on the lattice, essentially performing a discrete simulation 
of nonrelativistic quantum many-body theory. The resulting model falls in the class of quantum lattice gas 
automata, which were recently defined by MeyerQ in the context of the (1-1- l)-dimensional Dirac model. 
The exponential speedup-iof this alg»sithm on a quantum computer is a specific example of the general 
observation by Feynmanll3 and LloydlU that quantum mechanical systems can be simulated more efficiently 
on a quantum computer than on a classical computer. 

A quantum-computing device is composed of simple quantum elements such as particles with spin 1/2 
(q-bits) . The state space of the system at any fixed time is the tensor product of the Hilbert spaces of the 
states of the elementary computational elements. Thus, for example, a system with m qphits has a state 
space of dimension 2™. At each time step, some small number of q-bits (usually two or threeEj) are subjected 
to a unitary time evolution, described by acting with a unitary matrix on the Hilbert space of the affectc 
elements. Quantum computers have recently become of great interest because of the result due to Shorl 
that it is possible to factor large integers on a quantum computer in polynomial time, a procedure thought 
to be impossible on a classical computer. 

In order to implement the many-body simulation described in the beginning of this section on a quantum 
computer, it is necessary to make some restrictions on the behavior of the many-body wavefunction under 
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exchange of particles. The example system in Eq. (IV.5) describes two particles moving in one dimension 



without interacting. In this model, both particles can be moving in the same direction from the same lattice 
site at a given point in time. We can modify this model slightly to give the particles exclusionary (Fermi) 
statistics by making the transition matrix at x = y force the two particles to move in different directions. 
This corresponds to introducing a contact interaction between the two particles when they move within a 
single lattice distance. By making the initial conditions antisymmetric under exchange of x and y, we have 
a simulation of two nonrelativistic fermions moving in one dimension. Alternatively, we could symmetrize 
the wavefunction and we would have a simulation of "hard bosons" which obey Bose statistics but which 
cannot occupy the same lattice site. Either of these approaches naturally generalize to arbitrary numbers of 
particles and arbitrary dimensions. For the remainder of the discussion we assume that the particles, obey 
Bose statistics. The issue of implementing fermionie. systems on a quantum computer is more subtlellj, and 
has recently been addressed by Abrams and LloydtZl. 



10 



In previous sections we discussed the motion of a single particle, with a wave function ■0fc(x, t). Now, 
we would like to consider the state space for a quantum system of many particles. A natural basis for the 
Hilbert space of such a system is the set of states in the fermionic Fock space associated with the spatial 
lattice; such states are identified by a set of occupation numbers Sfc(x) (taking values or 1) for each possible 

particle position x and internal index k. The Hilbert space of the model is thus 2™' dimensional, where m 
is the number of possible values of the internal index, and is the number of lattice sites. For example, a 
basis vector of the state space for a one-dimensional system with 4 lattice sites x = 1,2,3,4 might be given 

by 

\s) - |(S2(1), si(l)), . . . , (S2(4), si(4))) = 1(0, 1), (0, 0), (0, 0), (1, 1)) (IV.6) 

where each ordered pair corresponds to the occupation numbers at a given lattice site. Thus, this state 
corresponds to the configuration where a single particle is at x = 1, with fc = 1, and both particle positions 
at a; = 4 are filled. 

The state of the quantum system at any given value of the discrete time parameter t is given by a vector 

s 

where the sum is taken over all basis vectors of the Hilbert space. This state is defined by the coefficients 
Cs{t). In the quantum computing paradigm, this corresponds to the state space of ml'^ independent q-bits. 

We will now define a quantum lattice gas automaton by defining a dynamics on the quantum state space. 
The dynamics of the quantum lattice-gas will be described in two steps, just as in classical lattice-gas 
automaton models. First there is a collision step in which the particles at each lattice site interact. Then 
there is an advection step, describing the propagation of the particles in the directions associated with the 
vectors c^. Each of these steps is described in the quantum system by a unitary transfer matrix acting on 
the state space of the system. The total dynamics can then be described by the equation 

m + At)) = A-K-m)) 

The advection step simply corresponds to a permutation matrix A on the basis vectors described above, 
where each bit is moved forward in the direction corresponding to the appropriate vector Cfc. For example. 



acting on the state in Eq. ([V.6), the result of applying the advection operator would be 

A\s) = 1(0,1), (0,1), (1,0), (0,0)) 

where we assume periodic boundary conditions on the lattice. The particle which was at lattice site x — I 
has been advected to x = 2,fc = 1, and the two particles which were at x = 4 have moved to a; = 3 and 
x = l. 

We now consider the collision part of the time development rule. The collision process is defined by a 
single unitary 2™ by 2™ matrix T, which acts separately on the quantum bits associated with each lattice 
site. Thus, the state of the system is transformed by the unitary matrix 

K ^ T(g)T(g) ■■■(g)T 

given by the I'^-iold tensor product of T. We would like the collision matrix T to have the property that it 
conserves particle number. Thus, this matrix is block diagonal in the subspaces of the Hilbert space at each 
lattice site corresponding to a fixed particle number. 

We have now defined a discrete model for quantum many-particle systems. To understand the behavior of 
this model in the continuum limit, let us consider the behavior when the number of particles in the system 
is relatively small compared to the number of lattice sites. In this case, at most lattice sites the number 
of particles present will be either or 1. This part of the dynamics, which describes the free propagation 
of single particles, is described by the part of T in the single-particle Hilbert space. However, because this 
is a unitary matrix, generically the dynamics described by this transition matrix is precisely that which 
we studied in the previous sections, and corresponds to a nonrelativistic particle propagating according to 
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the Schrodinger equation. Thus, for relatively sparse systems, this quantum lattice-gas model simiilatcs a 
system of many nonrelativistic particles whose free propagation is given by the Schrodinger equation. The 
remaining parts of the collision matrix T describe a contact interaction between the various particles. 

Let us now discuss the computational complexity of the quantum algorithm. To implement the advection 
transformation by using quantum computing elements, it is only necessary to perform a series of exchanges of 
the values of the quantum bits representing the particle occupation numbers. The number of such exchanges 
is essentially equal to the number of bits, ml'^ (recall that on a Euclidean lattice of dimension d, m = 2d, 
so that for example if = 3, the advection operation can be implemented in approximately 6/^ quantum 
operations). 

The matrix T acts on the Hilbert space associated with a subset of m of the q-bits in the system. Counting 
degrees of freedom, generically such a matrix can be implemented with approximately 2^™/15 elementary 
quantum operations on pairs of q-bits. For example, in a 3D system, it would take on the order of 300 
qiiantum operations to implement each T matrix, so that the number of computational steps needed to 
perform the transformation by K would be around 300/^. Note, however, that the part of T which acts on 
the multiple particle Hilbert space simply changes the phases of a delta function type interaction between 
the particles. Since these phases may not affect the results in most problems of physical interest, these 
components of T can be arbitrary. Thus, in practice we need only find a combination of operations on q-bits 
which will give a matrix T which preserves particle number and gives the desired symmetry properties and 
eigenvalues, reducing the number of steps needed significantly below 300. 

Combining these observations, we see that this model can be simulated with on the order of l'^ elementary 
quantum computations at each time step (or on the order of 1 if we are using a quantum computing system 
which allows parallel computation). Since this system automatically contains the multi-particle wave function 
for all possible particle numbers, we have achieved an exponential increase in speed over what was possible 
on a classical computer. 

The system as defined so far only includes interactions between the particles in the form of delta function 
interactions parameterized by the components of T in the multiple particle space. We can introduce an 
arbitrary interparticle potential V(x, y) by hand, by multiplying the wave function at each time step by the 
tensor product over all pairs of q-bits 

U = (S)J/i,7.X.V 



where the matrix 



Ui. 



/ 1 

10 

1 

\ e-'^'^(^'y) 



acts on the Hilbert space associated with the q-bits .Si(x) and Sj(y), changing the phase of the wavcfunction 
only in the component where both q-bits have the value 1. Implementing this interparticle potential will 
take on the order of m^Z^"^ quantum computations for each time step. Although this significantly increases 
the computational complexity of the quantum algorithm, this is still exponentially faster than the analogous 
classical algorithm, since the particle number n does not affect the complexity. Note that, unlike the rest of 
the algorithm, the implementation of interparticle potentials involves nonlocal interactions on the lattice. 

To clarify the discussion, we consider a simple example of a collision matrix T. For a many-body system 
in one dimension, at each lattice site the collision matrix T is a 4-by-4 matrix, acting on the Hilbert space 
with basis |(0,0)), |(0, 1)), |(1, 0)), |(1, 1)). Since we are assuming that particle number is conserved and that 
the dynamics is symmetric under left-right reflection, the matrix T is of the form 

/ Q 

a 6 

6 a 

V /3. 



where a,l3,a,b are complex numbers satisfying \a\' 



+ |6p = 1 and ab + ab = 0. By a simple 



global phase redefinition, we can choose a = 1. The part of T in the single-particle Hilbert space is precisely 
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the form of the coUision matrix S from Sec. From the eigenvalues of this matrix we can determine the 
mass of the free particles in the model. Finally, there is a single parameter (3 which describes the phase 
with which two particles "bounce" . In this simple one-dimensional model, there is therefore little freedom 
in choosing the particle interaction. In higher dimensions there would be nontrivial phases describing delta 
function interactions between up to 2d particles. 

One major concern in the implementation of any algorithm is the issue of precision. This problem is 
particularly acute on a quantum computer, where each quantum operation involves acting on the state 
with a unitary transformation which can only be controlled up to some finite precision. Furthermore, on 
a quantum computer there is the related but distinct problem of decoherence which must be addressed in 
order for any quantum computation to be feasible. There has been a great deal of workjxecently describing 
how these problems can be solved using dynamical quantum error correction methods£j. Without going 
into this issue in depth, we make the simple observation that even without error correction, if the precision 
of each quantum operation is 1 — l/t, then at each time step the error in the wave function will take a 
random step in the Hilbert space with size l/t. Only after on the order of operations will this error 
become significant. Thus, if we could achieve a precision better than 10~^, we could perform 10^° quantum 
operations successfully, which would allow us to simulate for example-an interacting 3D system on a lattice 
with size of order 20'^. With the error correction schemes described inll3, there is in principle no upper bound 
on how large a system could be simulated, other than the size of the quantum computer which could be built 
to perform the simulation. 

Finally, we consider the issue of measurement in quantum lattice gases. In a classical lattice gas, hydro- 
dynamic quantities, such as mass and momentum density, are obtained by averaging particles' mass and 
momentum over blocks in space and/or time. In a typical lattice-gas simulation, this is done from time to 
time to obtain the macroscopic variables of interest. The process of measuring these quantities is purely 
passive - that is, their measurement does not affect the subsequent dynamical evolution at all. In contrast, 
the analogous operation for a quantum lattice gas would involve occasionally measuring the state of some 
subset of the q-bits in the system, thus collapsing the quantum wavefunction onto the eigenstates of the 
(space and/or time) block number operator. The set of quantities which are accessible through this type of 
simulation are rather different from those accessible through simulation methods on a classical computer. 
For example, the dynamics of the system defines an effective Hamiltonian which is an approximation to the 
Hamiltonian of the many-body quantum system being simulated, however the spectrum of this Hamiltonian 
is not directly amenable to measurement. Instead, the types of observables which can be measured in the 
simulation are precisely equivalent to the types of observables which can be measured in an actual interacting 
quantum system. For example, a typical experiment might be to initialize the system in a particular known 
state at time t = 0, and to ask for the probability p at time t that there is a particle in a region of space dx^ . 
Just as in the physical quantum system, we can ask such a question of our simulation; we can perform the 
experiment a number of times, and each time we will find a particle with probability p. To actually compute 
p to some degree of accuracy requires repeating the experiment a number of times. 

V. NUMERICAL RESULTS 

To test the algorithm, we consider the dispersion relation of plane waves in periodic geometry in two 
dimensions. We consider a periodic grid with dimensions N hy N , and initialize it with a plane wave of the 
form 

ipj (x, 0) ~ ^ 6xp (ik • X — iujt) 

for j — 1, ... ,4, where 

k = 27r + lyf) , 

where Ix and ly are integers. Choosing units where the spatial dimensions are of unit length, we have 

e = and At ^ = 
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We evolve this initial condition in time, using Eq. ( |A.1| ) with ji = —i (hence m = 2) for the collisions. 
Every four time steps, we measure the inner product of the wave function with its initial condition, 



N2 

X 

The result should go like exp(—iujt), so the ratio of two successive values of this quantity is 
and hence the frequency is given by 

"^iA^^H S{t) )■ 

For a given wavevector k, we measure this frequency at many time steps t and take an average. 
We expect the evolution of the system to be governed by the Schrodinger equation, 

2m 



Since m — 2, this leads to the dispersion relation 



4 ■ 



We performed a series of simulations, where we considered wavenumbers = 31 and ly = /, where 

Z e {1, . . . , 12}. The points plotted in Fig. |l| show the measured frequency w as a function of |k| = 271^1"^ + l"^. 

The solid curve is |kp/4. It is evident that the agreement is excellent in the "hydrodynamic" limit of small 
|k|, but degrades due to lattice artifacts of order |k|Aa; at higher wavevector magnitudes. To demonstrate 
this, we include data for N — 256 (gray points) and N — 512 (black points). It is evident that the dispersion 
relation is valid for higher wavenumbers on the larger lattice. 
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FIG. 1. Plane- Wave Dispersion Relation is shown for = 256 (gray points) and = 512 (black points). 

VI. CONCLUSIONS 

We have considered a very general class of lattice models satisfying unitary time-evolution rules. We have 
shown that generic models of this type describe the evolution of a nonrelativistic particle according to the 
Schrodingcr equation in an arbitrary number of dimensions. These models can naturally be used to construct 
quantum lattice gases describing nonrelativistic many-body physics in an arbitrary number of dimensions. 
It is straightforward to include an arbitrary interparticle potential into these models. 

There are many ways in which this work could be extended. Numerical simulations could be performed in 
an arbitrary number of dimensions with multiple particles and with nontrivial spatially dependent potentials, 
and the results checked in analytically tractable cases. Further analysis is needed to understand the behavior 
of the system in the regime where particles are dense. It might also be interesting to consider more general 
collision rules which create and destroy particles, possibly including antiparticles with separate quantum 
numbers. 

Of course, the actual implementation of quantum lattice-gas models on quantum computing devices is 
something which may not be possible for many years, if ever. However, these lattice models give a simple 
framework with which to study problems in many-body theory. Furthermore, the methods described here 
arc also quite practical for simulating systems of a few particles on a classical computer. It may be that the 
exact unitarity of these models at a microscopic level will make them more stable and possibly more useful 
than currently used discrete methods such as finite-difference approaches. 
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APPENDIX A: SCHRODINGER EQUATION IN 2D 

We now present the formalism described above explicitly in two dimensions. The matrices D and X are 
given by 



D = 

















1 














1 














-1 
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X = 



(1/2 


Vl/2 



1/2 

1 

-1/2 



1/2 
]_ 


1/2 



1/2 


!_ 

-1/2/ 



This gives us the colhsion matrix 



S = X-'DX 



/ H + I t-i + 1 /I — 3 fi 



1 

^ + 1 /i + 1 /i+1 /i — 3 
^ — 3 /i + 1 /i+1 /i + 1 
V/i+l /i — 3 /i + 1 /i + 1 



(A.l) 



We can compute 



= XC^X-^ = 
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V2 
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V2 









/ ^ 



= XC^X-^ = 



and thus 
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1 




4^ -4^ 

1 

\/2 



1 
\/5 
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" (1-m)V2 




]_ 

2v/2 



1 



/ 



= XC^X-^ = 

Combining these matrices together we find 

\ 
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(l-p)\/2 
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2V2 
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2^2 
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' 2m 



dtC 
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00/ 



die 



/ 2m 
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0/ 
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0^ 
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1m 








2m 








Vo 








0^ 



dxdyC, 



with m described as in Eqs. ([II. 2) and (111.3), with d = 2. 

As predicted by the general discussion above, the total amplitude contained in the first component of C, 
satisfies a Schrodinger equation 

9,vE'(x,i) =1^(92 + a2)vE'(x,i) 



where 
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*(x,i) = ^[V'l(x,t)+V2(x,i)+V3(x,i)+V4(x,t)]. 

It is interesting to note that while the variation of the fourth component of ( contains an oscillating phase, 
and thus has no interesting behavior on the time scale of interest, the second and third components obey 
separate second-order differential equations analogous to the Schrodinger equation, but without rotational 
invariance. 



APPENDIX B: SCHRODINGER EQUATION IN 3D 



Using the above formalism in 3D, 



have 



X = 
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/Lt + 1 /X + 1 H+l /X + 1 + 1 /Lt — 5 
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From these matrices we find that the total amphtude 



«'(x, t)=^^—^ [V'i(x, t) + V2(x, t) + ^/-slx, t) + ^4(x, t) + V5(x, t) + ?A6(x, i)] . 



V6 

satisfies the Schrodinger equation 



a^M^lx, t) = + a,2 + dl)-^{^, t) 



where as usual m is related to /i through Eqs. (III. 2) and (III. 3), with d — Z. 
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